מעבדה 1 - עיבוד תמונה
Rom Hirsch and Yarom Swissa
import cv2 as cv
import sys
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import NoNorm
from PIL import Image
import numpy.matlib
import pandas as pd
import math
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from scipy import signal
from skimage import io
#%% functions
def tableDf(df):
fig, ax = plt.subplots()
# hide axes
fig.patch.set_visible(False)
ax.axis('off')
ax.axis('tight')
ax.table(cellText=df.values, colLabels=df.columns, loc='center')
fig.tight_layout()
def imageShow(im, title = "image"):
"""display image"""
cv.imshow(title, im)
cv.waitKey(0)
cv.destroyAllWindows()
def uint8c(x):
"""custom uint8 - (if val < 0 then val = 0)"""
x2 = np.copy(x)
x2[x2 < 0] = 0
x2[x2 > 255] = 255
return np.uint8(x2)
def plotPlusImage(axP, axI, axN, im, title):
axP.plot(im[1, :])
axP.set_title(title)
axI.imshow(im, cmap="gray", norm=NoNorm())
axI.set_title(title)
imgMin = np.min(im.flatten())
imgMax = np.max(im.flatten())
normalImg = np.dot(255, (im - imgMin) / (imgMax - imgMin))
axN.imshow(normalImg, cmap="gray")
axN.set_title(f"normalize {title}")
#axN.plot(normalImg[0, :])
# plot stem with color
def stem_plot(n, val, color):
markerline1, stemlines1, baseline1 = plt.stem(n, val)
plt.setp(markerline1, 'markerfacecolor', color)
plt.setp(stemlines1, linestyle="-", color=color, linewidth=2)
# Create Stem plot with color
def plotStem(title, ylabel, xlabel, color, x, y):
plt.title(title)
plt.ylabel(ylabel)
plt.xlabel(xlabel)
stem_plot(x, y, color)
#plt.show()
#padded zeros to array , params (array ,size of the array after the padded)
def padded_zeros(arr, N): # padded zeros
padded_array = np.zeros(N)
padded_array[:len(arr)] += arr
return padded_array
def add_sp(img, probability = 0.1):
noise = np.random.uniform(size = img.shape)
salt = noise > 1 - 0.5 * probability
pepper = noise < 0.5 * probability
return np.clip(img.astype('float64') + salt * 255 - pepper * 255, 0, 255).astype('uint8')
salt_pepper_noisy = add_sp(image)
def random_gaussian(N = 100000, mu = 0, sigma = 1):
x = np.zeros(N)
for i in range(0, N):
x[i] += (np.random.normal(0, 1) * sigma) + mu
return x
def conv2d(image, filter, mode = 'full', add_offset = 0, clip_enabled = False, clip_min = 0, clip_max = 255):
pad_amnt_y = len(filter) - 1
pad_amnt_x = len(filter[0]) - 1
padded_img = np.zeros((len(image) + 2 * pad_amnt_y, len(image[0]) + 2 * pad_amnt_x))
padded_img[pad_amnt_y : len(padded_img) - pad_amnt_y, pad_amnt_x : len(padded_img[0]) - pad_amnt_x] = image.astype('float64')
output = np.zeros((len(image) + pad_amnt_y * 2, len(image[0]) + pad_amnt_x * 2))
for x in range(0, len(output[0]) - len(filter[0]) + 1):
for y in range(0, len(output) - len(filter) + 1):
output[y][x] = np.sum(padded_img[y: y + len(filter), x : x + len(filter[0])] * filter)
output_float = output[:len(output) - pad_amnt_y, :len(output[0]) - pad_amnt_x]
output = output_float + add_offset
if (clip_enabled):
output = np.clip(output, clip_min, clip_max)
return output.astype('uint8'), output_float
def SquareInSquare(color_backgroud, color_object, res = False, range_background = 400):
black_with_square = np.ones([range_background, range_background]) * color_backgroud
black_with_square[150:250, 150:250] = color_object
black_with_square = np.uint8(black_with_square)
if res:
Cw = (color_object - color_backgroud) / color_backgroud
return black_with_square, Cw
else:
fig = plt.figure(figsize=(10, 10))
plt.subplot(211)
plt.imshow(black_with_square, cmap="gray", norm=NoNorm())
plt.title(f"SquareInSquare Bg = {color_backgroud} Bo = {color_object}")
plt.subplot(212)
vals = black_with_square.flatten()
plt.hist(vals, 256, [0, 256])
Cw = (color_object - color_backgroud) / color_backgroud
plt.title(f"Waber Contrast is:{Cw}")
#plt.show()
SquareInSquare(220,50)
num_stripes = np.arange(0, 256, int(256/8))
print(num_stripes)
[ 0 32 64 96 128 160 192 224]
def stripes_uint8(height, width, num_stripes):
if np.isscalar(num_stripes):
num_stripes = np.arange(0, 256, int(256/num_stripes))
width_stripe = math.floor(width/len(num_stripes))
vector_row = np.ones(width)
vector_row = np.array([vector_row[(i * width_stripe): ((i + 1) * width_stripe)] * v for i, v in enumerate(num_stripes)])
vector_row = vector_row.flatten()
return np.uint8(np.matlib.repmat(vector_row,height,1))
#1.
height = 200
width = 300
num_stripes = [0, 50, 100, 150, 200, 250]
plt.figure()
plt.imshow(stripes_uint8(height, width, num_stripes), cmap="gray")
plt.title("Stripes")
Text(0.5, 1.0, 'Stripes')
N = 256
Fs = N
Ts = 1/Fs
fx = 5
t = np.linspace(0, N-1, N) * Ts
sx = N
imX = np.zeros((N, N))
imX = np.matlib.repmat(50*np.sin(2*np.pi*fx*t), N, 1)
# plot
fig, axs = plt.subplots(3, 4,figsize=(16,13))
plotPlusImage(axs[0, 0], axs[1, 0], axs[2, 0], uint8c(np.abs(imX)), 'A line of uint(abs(imX))')
plotPlusImage(axs[0, 1], axs[1, 1], axs[2, 1], uint8c(imX), 'A line of uint(imX)')
plotPlusImage(axs[0, 2], axs[1, 2], axs[2, 2], np.float32(imX), 'A line of imX')
plotPlusImage(axs[0, 3], axs[1, 3], axs[2, 3], uint8c(imX-np.min(imX.flatten())), 'A line of uint(imX-min(imX.flatten()))')
plt.show()
# question 1 -4
A = 7
f = 10
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
#t = np.linspace(0, L, int(L/Ts))
y = np.cos(2*np.pi*f*t)
N = len(t)
plt.figure(figsize=(15, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")
fs = 12
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
plt.subplot(223)
plt.plot(t, y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.ylim(-1,1)
plt.subplot(224)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")
Text(0.5, 1.0, '(abs(fft(y)) fs = 12')
#%% 5
A = 7
f = 10
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
y = uint8c((1 + y) * 128)
N = len(t)
plt.figure(figsize=(10, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = unit8(sin(2*pi*f*t)) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")
fs = 12
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
y = uint8c((1 + y) * 128)
N = len(t)
plt.subplot(223)
plt.plot(t,y)
plt.title(f"y = unit8(sin(2*pi*f*t)) fs = {fs}")
plt.subplot(224)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")
Text(0.5, 1.0, '(abs(fft(y)) fs = 12')
רק האות בתדירות דגימה של 200 הרץ עומד בתנאי נייקוויסט שהוא דורש תדר כפול מתדר ההאות כלומר המינימום שלא יווצרו לנו התחזויות הוא אות תגימה של 20
בדגימה של 200 קיבלנו בספקטרות תדר 10 כפי שאות המקור לאומת זאת בדגימה של 12 קיבלנו אות בתדר 2 זה נובע בגלל שהפרש בין התדגימה לאות הוא 2 לכן ציפינו לקבל אות בתדירות 2
בשאלות 5 כפי שציפינו נוסף רכיב דיסי לספקטרום בגלל הנירמול בין 0-255
f = np.array([1, 2, 3, 1, 2])
h = np.array([1, 1])
g = np.convolve(f, h)
n = np.arange(0, 4, 1)
m1 = len(f)
m2 = len(h)
f = padded_zeros(f, m2+m1-1)
h = padded_zeros(h, m2+m1-1)
plt.figure(figsize=(10, 10))
plt.subplot(311)
plotStem("f[n]", "amplitude", "n", 'red', range(m1+m2-1), f)
plt.subplot(312)
plotStem("h[n]", "amplitude", "n", 'green', range(m1+m2-1), h)
convXY = np.convolve(f, h)
plt.subplot(313)
plotStem("conv(f*h)", "amplitude", "n", 'blue', range(convXY.shape[0]), convXY)
plt.show()
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.
# 1.4 - q6
f = np.array([1, 2, 3, 1, 2])
h2 = np.zeros(f.shape[0]*2)
m1 = len(f)
m2 = len(h2)
h2[[0,f.shape[0]]] = 1
f = padded_zeros(f, m2)
plt.figure(figsize=(10, 10))
plt.subplot(311)
plotStem("f[n]", "amplitude", "n", 'red', range(m2), f)
plt.subplot(312)
plotStem("h[n]", "amplitude", "n", 'green', range(m2), h2)
g2 = np.convolve(f,h2)
plt.subplot(313)
plotStem("conv(f*h)", "amplitude", "n", 'blue', range(m2), g2[0:m2])
plt.show()
/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True. /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:41: UserWarning: In Matplotlib 3.3 individual lines on a stem plot will be added as a LineCollection instead of individual lines. This significantly improves the performance of a stem plot. To remove this warning and switch to the new behaviour, set the "use_line_collection" keyword argument to True.
# Gaussian Noise
mu = np.array([0, 3, 5, 10])
sigma = np.array([0.5, 2.5, 9, 10])
noise = np.array([np.random.normal(mu[i], sigma[i], 1000) for i in range(mu.shape[0])])
plt.figure(figsize=(14, 14))
for i in range(noise.shape[0]):
plt.subplot(4, 1, i+1)
plt.plot(noise[i])
plt.title(f"noise {i}")
plt.xlabel("samples", horizontalalignment='right', x=1.0)
plt.ylabel("amp")
plt.show()
plt.figure(figsize=(15, 15))
for i in range(noise.shape[0]):
plt.subplot(4, 1, i+1)
counts, bins = np.histogram(noise[i])
plt.hist(bins[:-1], bins, weights=counts)
plt.title(f"noise {i}")
plt.xlabel("samples", horizontalalignment='right', x=1.0)
plt.ylabel("amp")
plt.show()
8 הוספת רעש גאוסי והורדת הרעש ע"י מסנן
#%% 1.4.4
h1 = np.array([0.2]*5)
A = 7
f = 10
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
y = np.cos(2*np.pi*f*t)
n1 = np.random.normal(mu[0], sigma[0], len(y))
f1 = n1 + y
g1 = np.convolve(f1,h1)
plt.figure(figsize=(20,15))
plt.subplot(413)
plt.plot(f1)
plt.title("f1[n]")
plt.subplot(411)
plt.plot(n1)
plt.title("n1[n]")
plt.subplot(412)
plt.plot(y)
plt.title("y[n]")
plt.subplot(414)
plt.plot(g1)
plt.title("g1[n]")
Text(0.5, 1.0, 'g1[n]')
מסקנות : הרעש שהתווסף הוא רעש גאוסי אשר מתפרס על כל התדרים ולכן לאחד שהעברנו את האות במסנן מעביר נמוכים נפתרנו מרעש בתדרים הגבוהים אך עדיין נשאר הרעש בתדרים הנמוכים
ממדי האות לאחר הקונבולוציה יצאו הגודל של אות המקורי ועוד הגודל של המסנן פחות אחד וזה מתאים לתאוריה
plt.figure(figsize=(20,15))
plt.subplot(311)
plt.title("h1[n]")
plt.plot(np.linspace(-fs/2, fs/2,h1.shape[0]),np.fft.fftshift(np.fft.fft(h1)))
plt.subplot(312)
plt.title("f1[n]")
plt.plot(np.linspace(-fs/2, fs/2,f1.shape[0]),np.fft.fftshift(np.fft.fft(f1)))
plt.subplot(313)
plt.title("g1[n]")
plt.plot(np.linspace(-fs/2, fs/2,g1.shape[0]),np.fft.fftshift(np.fft.fft(g1)))
/usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1317: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float) /usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1317: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float) /usr/local/lib/python3.7/dist-packages/matplotlib/cbook/__init__.py:1317: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)
[<matplotlib.lines.Line2D at 0x7f7a4ddd2ad0>]
קוד -המששמש לאפליקציה אינטרקטיבית לקביעת ניגודיות וובר
interactive application :
control keys:
object brightness: z - increase, x - decrease ;
background brightness: a - increase, s - decrease, g - display histogram + save state to the table q - exit
#%% interactive application : key: object z - increase, x - decrease ;
# background a - increase, s - decrease,
# g - display histogram + save state to the table, q - exit
bg=1
bo=1
d = {"brightness (bg)" : ["brightness (bo)", "JND(Cw)"]}
c = 1
# to
"""
while (True):
image,Cw = SquareInSquare(bg,bo,True)
cv.imshow(f"bg = {bg}, bo = {bo}",image)
k = cv.waitKey(0)
if k == ord('a'):
bg += 1
elif k == ord('s'):
bg = bg - 1
elif k == ord('z'):
bo += 1
elif k == ord('x'):
bo = bo - 1
elif k == ord('q'):
cv.destroyAllWindows()
break
elif k == ord('g'):
#addToTable()
if f"{bg}" in d.keys():
d[f"{bg}({c})"] = [bo,Cw]
else:
d[f"{bg}"] = [bo,Cw]
SquareInSquare(bg,bo)
else:
continue
cv.destroyAllWindows()
df = pd.DataFrame(data=d)
print(df.to_markdown()) # require pip install tabulate
"""
SquareInSquare(220,180)
התוצאות שלנו המניסוי: במעבדה בכיתה באור דלוק
#%%resulte of the test:
d = {"brightness (bg)" : ["brightness (bo)", "JND(Cw)"], '40': [39, -0.025],'40(1)': [41, 0.025],'140': [140, 0.0],'140(1)': [141, 0.007142857142857143],'220': [221, 0.004545454545454545],'220(1)': [220, 0.0]}
df = pd.DataFrame(data=d, index=['brightness(bo)', 'Cw'])
print(df.to_markdown(index=False)) # require pip install tabulate
| brightness (bg) | 40 | 40(1) | 140 | 140(1) | 220 | 220(1) | |:------------------|-------:|--------:|------:|-------------:|-------------:|---------:| | brightness (bo) | 39 | 41 | 140 | 141 | 221 | 220 | | JND(Cw) | -0.025 | 0.025 | 0 | 0.00714286 | 0.00454545 | 0 |
תוצאות מבדיקה נוספת בבית בתנאים אחרים:
d = {'brightness (bg)': ['brightness (bo)', 'JND(Cw)'],
'40': [41, 0.025],
'140': [139, -0.007142857142857143],
'140(1)': [142, 0.014285714285714285],
'220': [219, -0.004545454545454545],
'220(1)': [221, 0.004545454545454545]}
df = pd.DataFrame(data=d, index=['brightness(bo)', 'Cw'])
print(df.to_markdown(index=False)) # require pip install tabulate
| brightness (bg) | 40 | 140 | 140(1) | 220 | 220(1) | |:------------------|-------:|-------------:|------------:|-------------:|-------------:| | brightness (bo) | 41 | 139 | 142 | 219 | 221 | | JND(Cw) | 0.025 | -0.00714286 | 0.0142857 | -0.00454545 | 0.00454545 |
גרף תאורתי של רגישות :
התדר הדומיננטי בגל ריבועי הוא תדר מאוד נמוך לכן הC
The dominant frequency of square wave is near DC, so by looking at the theoretical graph we would expect CW around 30 A
כושר ההבחנה תלוי גם בתאורה בחדר בסוג המסך ובצופה לכן יתקבל תוצאות שונות לכל בדיקה במידה ותנאים אלו ישתנו
f, ax = plt.subplots(4, 3,figsize=(25,25))
img=stripes_uint8(height, width, 8)
ax[0, 0].imshow(img,cmap="gray")
ax[0, 0].set_title("8 bands Image")
ax[0, 0].axis('off')
print(img.flatten()[-1])
ax[0, 1].hist(img.flatten(), 256, [0, 256])
ax[0, 1].set_title("8 bands Histogram")
ax[0, 2].plot(img[0,:])
ax[0, 2].set_title("the 1st line")
img=stripes_uint8(height,width,16)
ax[1, 0].imshow(img,cmap="gray")
ax[1, 0].set_title("16 bands Image")
ax[1, 0].axis('off')
ax[1, 1].hist(img.flatten(), 256, [0, 256])
ax[1, 1].set_title("8 bands Histogram")
ax[1, 2].plot(img[0,:])
ax[1, 2].set_title("the 1st line")
img=stripes_uint8(height, width, 32)
ax[2, 0].imshow(img,cmap="gray")
ax[2, 0].set_title("32 bands Image")
ax[2, 0].axis('off')
ax[2, 1].hist(img.flatten(), 256, [0, 256])
ax[2, 1].set_title("32 bands Histogram")
ax[2, 2].plot(img[0,:])
ax[2, 2].set_title("the 1st line")
img=stripes_uint8(height, width, 256)
ax[3, 0].imshow(img,cmap="gray")
ax[3, 0].set_title("Ramp Image")
ax[3, 0].axis('off')
ax[3, 1].hist(img.flatten(), 256, [0, 256])
ax[3, 1].set_title("Ramp Histogram")
ax[3, 2].plot(img[0,:])
ax[3, 2].set_title("the 1st line")
plt.show()
img=stripes_uint8(height, width, 8)
noise = np.uint8(np.matlib.repmat(np.random.normal(0, 20, img.shape[0]),img.shape[1],1))
224
הוספת רעש:
img=stripes_uint8(height, width, 8)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)
f, ax = plt.subplots(4, 3,figsize=(25,25))
img=stripes_uint8(height, width, 8)
img = img + gauss
ax[0, 0].imshow(img,cmap="gray")
ax[0, 0].set_title("8 bands Image")
ax[0, 0].axis('off')
print(img.flatten()[-1])
ax[0, 1].hist(img.flatten(), 256, [0, 256])
ax[0, 1].set_title("8 bands Histogram")
ax[0, 2].plot(img[0,:] )
ax[0, 2].set_title("the 1st line")
img=stripes_uint8(height,width,16)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)
img = img + gauss
ax[1, 0].imshow(img,cmap="gray")
ax[1, 0].set_title("16 bands Image")
ax[1, 0].axis('off')
ax[1, 1].hist(img.flatten(), 256, [0, 256])
ax[1, 1].set_title("8 bands Histogram")
ax[1, 2].plot(img[0,:])
ax[1, 2].set_title("the 1st line")
img=stripes_uint8(height, width, 32)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)
img = img + gauss
ax[2, 0].imshow(img,cmap="gray")
ax[2, 0].set_title("32 bands Image")
ax[2, 0].axis('off')
ax[2, 1].hist(img.flatten(), 256, [0, 256])
ax[2, 1].set_title("32 bands Histogram")
ax[2, 2].plot(img[0,:])
ax[2, 2].set_title("the 1st line")
img=stripes_uint8(height, width, 256)
row,col= img.shape
mean = 0
sigma = 20
gauss = np.random.normal(mean,sigma,(row,col))
gauss = gauss.reshape(row,col)
img = img + gauss
ax[3, 0].imshow(img,cmap="gray")
ax[3, 0].set_title("Ramp Image")
ax[3, 0].axis('off')
ax[3, 1].hist(img.flatten(), 256, [0, 256])
ax[3, 1].set_title("Ramp Histogram")
ax[3, 2].plot(img[0,:])
ax[3, 2].set_title("the 1st line")
plt.show()
img=stripes_uint8(height, width, 8)
noise = np.uint8(np.matlib.repmat(np.random.normal(0, 20, img.shape[0]),img.shape[1],1))
266.1002506128287
class d2signal:
def __init__(self, sx, sy, fs, fx, fy, amp=1, phase=0, func=np.sin):
self.n = N
self.amp = amp
self.fx = fx
self.fy = fy
self.fs = fs
self.func = func
self.sx = np.array(np.linspace(0, sx-1, sx)).reshape((1,-1))
self.sy = np.array(np.linspace(0, sy-1, sy).reshape((-1,1)))
self.img = self._genSignal()
def _genSignal(self):
Ts = 1/self.fs
wx = 2 * np.pi * self.fx * Ts;
wy = 2 * np.pi * self.fy * Ts;
self.img = 127 - self.amp * self.func(wy * self.sy + wx * self.sx)
self._signalUint8()
return self.img
def _signalUint8(self):
self.imgUint8 = uint8c(self.img)
def imshow(self,title="image", Uint8 = False):
if Uint8:
img = self.imgUint8
else:
img = self.img
plt.figure(figsize=(15,15))
plt.title(title)
plt.imshow(img, cmap="gray",norm=NoNorm())
plt.show()
def normalize(self):
imgMin = np.min(self.img.flatten())
imgMax = np.max(self.img.flatten())
self.old = self.img.copy()
self.img = np.dot(255, (self.img - imgMin) / (imgMax - imgMin))
self._signalUint8()
def fft(self, title="fft", Uint8 = False):
plt.figure(figsize=(15,15))
if Uint8:
img = self.imgUint8
else:
img = self.img
plt.title(title)
plt.imshow(np.fft.fftshift(abs(np.fft.fft2(img))), cmap="gray")
plt.show()
def imshow3d(self, Uint8 = False):
from matplotlib import cm
from matplotlib.ticker import LinearLocator
if Uint8:
img = self.imgUint8
else:
img = self.img
xx, yy = np.mgrid[0:img.shape[0], 0:img.shape[1]] # get an x y thing of the image
fig = plt.figure(figsize=(15,15))
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xx, yy, img, cmap=cm.gray,
linewidth=0, antialiased=False)
ax.set_zlim(0, 255)
ax.zaxis.set_major_locator(LinearLocator(10))
ax.zaxis.set_major_formatter(FormatStrFormatter('%d'))
fig.tight_layout()
plt.show()
I1 = d2signal(256,256,256,5,0,50)
I1.imshow("Amp = 50 , fx = 5",True)
#%%
I2 = d2signal(256,256,256,0,10,100)
I2.imshow("Amp = 100 fy= 10",True)
#%%
I3 = d2signal(256,256,256,5,10,127)
I3.imshow("Amp = 127 fy= 10 fx = 5",True)
מה שניתן לראות בתמונה 3 זה שיש לנו תדר שהוא בשני הצירים. לפי הספקטרום התדר ניתן לקבל את התדר המרחבי ואת הרכיבים בכל ציר
I1.fft()
I2.fft()
I3.fft()
I1.normalize()
I2.normalize()
I3.normalize()
I1.imshow("Amp = 50 , fx = 5",True)
I2.imshow("Amp = 100 fy= 10",True)
I3.imshow("Amp = 127 fy= 10 fx = 5",True)
I4 = d2signal(256,256,256,5,0,10)
I4.img = I1.img + I2.img + I3.img
I4.normalize()
I4.imshow("Amp = 50 , fx = 5",True)
I4.fft()
I4.imshow3d()
min_freq = 0.1 # cycles per image
N=400
min_amplitude = 1
max_amplitude = 127
max_freq = 40
f = np.exp(np.linspace(np.log(min_freq), np.log(max_freq), N).reshape((1,-1)))
a = np.exp(np.linspace(np.log(min_amplitude), np.log(max_amplitude), N).reshape((-1,1)))
img = 128 + a * np.sin(2 * np.pi * f)
plt.figure(figsize=(15,15))
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
<matplotlib.image.AxesImage at 0x7f7a4fbb1290>
filters + mse + snr + psnr:
url = r"https://placekitten.com/500/300"
image = io.imread(url,0)
image = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
#create filters
lpf = np.array([[1,1,1],[1,1,1],[1,1,1]]) / 9.0
hpf = np.array([[-1,-1,-1],[-1,8,-1],[-1,-1,-1]]) / 9.0
sharp_f = np.array([[-1,-1,-1],[-1,17,-1],[-1,-1,-1]]) / 9.0
laplacian_f = np.array([[0,-1,0],[-1,4,-1],[0,-1,0]])
im_lp, _ = conv2d(image, lpf)
im_laplacian, im_laplacian_float = conv2d(image, laplacian_f, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
im_hp, im_hp_float = conv2d(image, hpf, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
median_filt = cv.medianBlur(image, 3)
salt_pepper_noisy = add_sp(image, 0.1)
gaussian_noise = random_gaussian(N = (len(image) * len(image[0])), sigma = 25)
gaussian_noise = np.reshape(gaussian_noise, image.shape)
gaussian_noisy = np.clip((image.astype('float64') + gaussian_noise), 0, 255).astype('uint8')
salt_pepper_lpf, _ = conv2d(salt_pepper_noisy, lpf)
gaussian_lpf, _ = conv2d(gaussian_noisy, lpf)
salt_pepper_hpf, salt_pepper_hpf_float = conv2d(salt_pepper_noisy, hpf, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
gaussian_hpf, gaussian_hpf_float = conv2d(gaussian_noisy, hpf, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
gaussian_sharp, _ = conv2d(gaussian_noisy, sharp_f, clip_enabled = True)
salt_pepper_sharp, _ = conv2d(salt_pepper_noisy, sharp_f, clip_enabled = True)
salt_pepper_laplacian, salt_pepper_laplacian_float = conv2d(salt_pepper_noisy, laplacian_f, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
gaussian_laplacian, gaussian_laplacian_float = conv2d(gaussian_noisy, laplacian_f, add_offset = 127, clip_enabled = True, clip_min = 0, clip_max = 255)
salt_pepper_median = cv.medianBlur(salt_pepper_noisy, 3)
im_sharp, _ = conv2d(image, sharp_f, clip_enabled = True, clip_min = 0, clip_max = 255)
gaussian_median = cv.medianBlur(gaussian_noisy, 3)
salt_pepper_fixed, _ = conv2d(salt_pepper_median, sharp_f, clip_enabled = True)
# calc mse + snr + snr
img_ref = image
reference_max = np.max(img_ref).astype('float64') #
reference_avg = np.average(img_ref).astype('float64')
reference_variance = np.sum((image.astype('float64') - reference_avg) * (image.astype('float64') - reference_avg) / (len(image) * len(image[0])))
plt.figure(figsize=(15,15))
figure, axes = plt.subplots(3, 3,figsize=(15,15))
axes[0][0].imshow(img_ref, cmap = 'gray', vmin = 0, vmax = 255)
axes[0][0].set(title = 'original')
axes_index = 1
imgs = [im_lp, im_hp_float, im_sharp, im_laplacian_float, gaussian_noisy, gaussian_lpf, salt_pepper_noisy, salt_pepper_fixed]
labels = ['orig - lpf', 'orig - hpf', 'orig - sharp', 'orig - laplacian', 'gaussian noise', 'gaussian - lpf', 'sp - noise', 'sp - nmedian&nsharp']
snr_log = []
psnr_log = []
mse = []
for img in imgs:
# because the convolution increased the size of the image, ignore the borders
if ((len(img) != len(img_ref)) and (len(img[0]) != len(img_ref[0]))):
img_ = img[1:len(img) - 1,1:len(img[0]) - 1]
else:
img_ = img
noise = img_ref.astype('float64') - img_.astype('float64')
noise_avg = np.average(noise)
noise_variance = np.sum((noise - noise_avg) * (noise - noise_avg)) / (len(noise) * len(noise[0]))
snr_log.append(10*np.log10(reference_variance / noise_variance))
current_mse = np.sum(noise * noise) / (len(noise) * len(noise[0]))
mse.append(current_mse)
psnr_log.append(10*np.log10(reference_max * reference_max / current_mse))
axes[int(axes_index / 3)][int(axes_index % 3)].imshow(img_, cmap = 'gray', vmin = 0, vmax = 255)
axes[int(axes_index / 3)][int(axes_index % 3)].set(title = labels[axes_index - 1])
axes_index += 1
figure.tight_layout()
<Figure size 1080x1080 with 0 Axes>
כאשר מופיע רעש גאוסייני סדר הפעולות הסינון היעיל ביותר יהיה מסנן מעביר נמוכים ואז מסנן sharp
כאשר מופיע רעש מלח פלפל היעיל ביותר יהיה פילטר חציון
plt.figure(figsize=(15,10))
plt.xticks(np.arange(len(labels)),labels)
plt.plot(snr_log)
plt.plot(psnr_log)
plt.title('SNR, PSNR, log10 scale')
plt.legend(['SNR','PSNR'])
plt.show()
plt.figure(figsize=(15,10))
plt.xticks(np.arange(len(labels)),labels)
plt.plot(mse)
plt.title("MSE")
plt.show()
מספר הרמות המצג בתרשים פסי מאך ללא רעש מוגדר כך שככל שמספר רמות האפור גדול יותר העין מתפקדת כאינטגרטור ולכן יהיה יותר קשה להבחין בין פס לפס
לעומת זאת ככל שמספר רמות האפור קטן יותר העין מתפקדת הגוזר וניתן לראות את המעברים בין פס לפס בקלות יותר. העין יכולה להבחין בערך 100 רמות אפור
A = 1
f = 66.67
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
#t = np.linspace(0, L, int(L/Ts))
y = A * np.cos(2*np.pi*f*t)
N = len(t)
print(f"min: {np.min(y)} max: {np.max(y)}")
plt.figure(figsize=(15, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")
min: -0.517668669793734 max: 1.0
Text(0.5, 1.0, '(abs(fft(y)) fs = 200')
A = 1
f = 65.93
L = 1
fs = 200
Ts = 1/fs
t = np.arange(0,L-Ts,Ts)
#t = np.linspace(0, L, int(L/Ts))
y = A * np.cos(2*np.pi*f*t)
N = len(t)
print(f"min: {np.min(y)} max: {np.max(y)}")
plt.figure(figsize=(15, 10))
plt.subplot(221)
plt.plot(t,y)
plt.title(f"y = sin(2*pi*f*t) fs = {fs}")
plt.subplot(222)
plt.plot(np.linspace(-(fs/2), fs/2-1, y.shape[0]), np.fft.fftshift(abs(np.fft.fft(y))))
plt.title(f"(abs(fft(y)) fs = {fs}")
min: -0.9998507259473712 max: 1.0
Text(0.5, 1.0, '(abs(fft(y)) fs = 200')
סעיף ה: מהתוצאות ניתן לראות שכאשר דוגמים סינוס לא בכפולות שלמות, אנחנו מקבלים אות מזוגזג. כלומר, אנו דוגמים את הסינוס באפליטודה שונה בכל דגימה. כדי לקבל סינוס חלק, נדגום את האות בכפולה שלומה ובכך נדגום בכל מחזור את אותה אפליטודה של האות.
סעיף ב: הערכים שמפריעים לנו הם ערכים הגבוהים מ255, אם נוריד ערך קבוע מהמטריצה (שווה ערך לערך DC קבוע כלשהו) נוכל לטפל בבעיה ובכך למנוע ממצב בו מגיעים לערכי רוויה.
סעיף ג: הנירמול יגדיל את הטווח הדינמי ובכך לא יתרום מאחר והערכים נמצאים כבר ברוויה, במקרה הזה הטווח הדינמי מנוצל מאחר וקיימים ערכים מ0 ועד 255.
f = np.array([1, 2, 3, 4, 5])
print(f"f = {f}")
size = f.shape[0]
zerosPadding = 3
size = size + zerosPadding
h = np.array([1]+[0]*size)
numberDuplicate = 4
h = np.matlib.repmat(h,1,5).flatten()
m1 = len(f)
m2 = h.shape[0]
f = padded_zeros(f, m2+m1-1)
g = np.convolve(h,f)
print(f"h = {h}")
print("Result conv(f,h):")
print(g[:m2])
f = [1 2 3 4 5] h = [1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0] Result conv(f,h): [1. 2. 3. 4. 5. 0. 0. 0. 0. 1. 2. 3. 4. 5. 0. 0. 0. 0. 1. 2. 3. 4. 5. 0. 0. 0. 0. 1. 2. 3. 4. 5. 0. 0. 0. 0. 1. 2. 3. 4. 5. 0. 0. 0. 0.]
#%%
## create a vector of exponential values for frequencies
min_freq = 0.1 # cycles per image
N=400
min_amplitude = 1
max_amplitude = 255
max_freq = 50 # cycles per image
f = np.exp(np.linspace(np.log(min_freq), np.log(max_freq), N).reshape((1,-1)))
a = np.exp(np.linspace(np.log(min_amplitude), np.log(max_amplitude), N).reshape((-1,1)))
img = 128 + a * np.sin(2 * np.pi * f)
plt.figure(figsize=(15,15))
plt.subplot(131)
plt.title("Campbell-Robson")
plt.imshow(img, cmap='gray', vmin=0, vmax=255)
plt.subplot(132)
plt.title("Magnitude")
amp_img = np.uint8(np.matlib.repmat(a,1,img.shape[1]))
plt.imshow(amp_img, cmap='gray', vmin=0, vmax=255)
plt.subplot(133)
cos_img = uint8c(np.matlib.repmat( 255*np.cos(2 * np.pi * f)+128,img.shape[1],1))
plt.title("Cosine")
plt.imshow(cos_img, cmap='gray', vmin=0, vmax=255)
plt.figure(figsize=(10,10))
plt.subplot(211)
plt.title("Amplitude as a function of row nember")
plt.xlabel("Row")
plt.ylabel("Amp")
plt.plot(a)
plt.subplot(212)
plt.title("Frequency as a function of row nember")
plt.xlabel("Row")
plt.ylabel("Frequency")
plt.plot(f.T)
[<matplotlib.lines.Line2D at 0x7f7a4fed7590>]
fmax = 50hz fmin = 1hz AmpMax = 255 AmpMin = 1
a. LPF is used to remove high frequency noise and smooth pictures HPF is used for edge detection Median filter is most commonly used to remove salt and pepper noise b. filters are used in cascade when we want diffrent effects on the picture, for example, if we have and picture blurred filter with salt and pepper noise we would like to use median and HPF\sharpening filters in cascade c. for 2 cascade filters(MxM) on a picture of NxN the complexity would be O(M^4*N^4), we can use the properties of convolution and compute the convolution of the smaller matrices first. if N < M than we could compute the filters together and than use the product on the picture. d. a median filter cannot be used for the solution becuase it is not a filter that applied by convolution.
4.1 - 4.5 בוצע במהלך הדוח
5.1 - 5.5 בוצע במהלך הדוח
imfilter
הפונקציה אשר ניתן להשתמש בה באופן רב ממדי לעומת conv2 אשר רצה על דו מימד
1)From the warmout excercises we concluded that a transition between two distant colors is translated to a high frequency it is easy for our eye to notice. 2) Images with diffrent noises requires specific filters, if the wrong filter is used it can make the image worse. 3) when using filters in cascade it is important to notice in what order we use the filters, for example, for a blurred image with salt and pepper noise we must use median filter first and sharpening filter after, if we switch the order the image will be substantially worse. 4) using one byte per pixel can sometimes cause problems, for example when adding images that contains sines over dc of 128 we reach values bigger than 255, converting to double and normalizing back to 0-255 can solve the problem.